home *** CD-ROM | disk | FTP | other *** search
/ Amiga Collections: Taifun / Taifun 038 (1987-11-15)(Ossowski, Stefan)(DE)(PD).zip / Taifun 038 (1987-11-15)(Ossowski, Stefan)(DE)(PD).adf / BasicProgs / LeastSquare < prev    next >
Text File  |  1989-01-18  |  13KB  |  526 lines
  1. DIM P(20), r(20), v(20), X(20), Y(20), z(20,20)
  2. DEFDBL a-z
  3. DEFSNG PrintY
  4. SCREEN 2,640,200,3,2
  5. WINDOW 2,,,0,2
  6. WIDTH 80
  7. Round$ = "################.####"
  8. MenuFlag = 0 'for degree count
  9.  
  10. LOCATE 10,28
  11. PRINT  "Want instructions (Y/N)?"
  12. Instruct:
  13.   Want$ = INKEY$
  14.   IF Want$ = "" THEN Instruct
  15.   Want$ = UCASE$(Want$)
  16.   IF Want$ <> "Y" AND Want$ <> "N" THEN Instruct
  17.   CLS
  18.   IF Want$ = "N" THEN GOTO Start
  19.   
  20. IntroFlag = 0
  21. RESTORE
  22. PrintIntro:
  23.   FOR j% = 1 TO 23
  24.     READ d$
  25.    Over = INT((80 - LEN(d$))/2)
  26.     PRINT  SPACE$(Over) d$
  27.   NEXT
  28.   EndIntro:
  29.     IF INKEY$ = "" THEN EndIntro
  30.     IF IntroFlag = 0 THEN
  31.       CLS
  32.       IntroFlag = 1
  33.       GOTO PrintIntro
  34.     END IF
  35.    
  36.  
  37. Start:
  38.   CLS
  39.   WIDTH 70
  40.   NumPairs = 0
  41.   HowManyPlaces:
  42.   LOCATE 5,11
  43.   PRINT  "How many decimal places do you want to show - 2 is mimimum"
  44.   LOCATE 6,11
  45.   INPUT "(Enter your choice then press RETURN)    "; Places
  46.     IF Places < 2 THEN
  47.       PRINT  "  Everybody wants to be a comedian.  Try again."
  48.       PRINT 
  49.       GOTO HowManyPlaces
  50.     END IF
  51.     IF Places > 16 THEN
  52.       PRINT  "16 is as high as the machine goes.  Try again."
  53.       GOTO Start
  54.     END IF
  55.   Place$ = "." + LEFT$(Round$,Places)
  56.   Round$ = Round$ + Place$
  57.   Places = Places + 1
  58.     
  59. InputLoop:
  60.   PRINT 
  61.   NumPairs = NumPairs+1
  62.   COLOR 3,0
  63.     PRINT "    Pair"; NumPairs;
  64.   COLOR 4,0
  65.     PRINT  "  (C to start calculating - D to delete last pair)"
  66.   COLOR 1,0
  67.     PRINT "    Input the first number of the pair then press RETURN"
  68.   GOSUB GetIt
  69.   IF Final$ = "C" AND NumPairs > 2 THEN NumPairs = NumPairs-1: GOTO SetDegree
  70.   IF Final$ = "D" AND NumPairs > 1 THEN NumPairs = NumPairs-2: GOTO InputLoop
  71.   X(NumPairs) = Final
  72.   PRINT 
  73.   PRINT "    Input the second Number of the pair then press RETURN"
  74.   GOSUB GetIt
  75.   IF Final$ = "C" AND NumPairs > 2 THEN NumPairs = NumPairs-1: GOTO SetDegree
  76.   IF Final$ = "D" AND NumPairs > 1 THEN NumPairs = NumPairs-2: GOTO InputLoop
  77.   Y(NumPairs) = Final
  78.   GOTO InputLoop
  79.  
  80. SetDegree:
  81.   CLS
  82.   FOR j% = 1 TO NumPairs  'sort in ascending order
  83.     FOR k% = j% TO NumPairs
  84.       IF X(j%) > X(k%) THEN
  85.         SWAP X(j%), X(k%)
  86.         SWAP Y(j%), Y(k%)
  87.       END IF
  88.     NEXT
  89.   NEXT
  90.   Degree = NumPairs - 1
  91.   IF Degree > 9 THEN Degree = 9
  92.   FOR WhatDeg = 1 TO Degree:  GOSUB Calculate:  NEXT
  93.  
  94. PrintMenu:
  95.   MenuFlag = 1
  96.   COLOR 5,0
  97.   PRINT 
  98.   PRINT 
  99.   PRINT SPACE$(38)"MENU"
  100.   PRINT 
  101.   PRINT SPACE$(20)"1 - Drop this equation and find another"
  102.   PRINT SPACE$(20)"2 - Determine Y given X"
  103.   PRINT SPACE$(20)"3 - List all pairs"
  104.   PRINT SPACE$(20)"4 - Quit"
  105.   PRINT SPACE$(20)"5 - Graph"
  106.   COLOR 1,0
  107.   PRINT "                                                 
  108.   GetQ:
  109.     q$ = INKEY$
  110.     q = VAL(q$)
  111.     IF q<1 OR q>5 THEN GetQ
  112.   ON q GOTO Start, FindY, ListPairs, EndIt, GraphIt
  113.   GOTO PrintMenu
  114.  
  115. FindY:
  116.   WhatDeg = 0: INPUT "What degree equation do you want to use"; WhatDeg
  117.   GOSUB Calculate
  118.   NextY:
  119.     PRINT 
  120.     PRINT "(R to return to menu)"
  121.     INPUT "x = "; X$
  122.     X = VAL(X$)
  123.     X$ = UCASE$(X$)
  124.     IF X$ = "R" THEN GOTO PrintMenu
  125.     YValue = 0
  126.     FOR k% = 1 TO WhatDeg + 1
  127.       YValue = YValue + v(k%) * X ^ (k%-1)
  128.     NEXT
  129.     PRINT "y ="; YValue
  130.     GOTO NextY
  131.  
  132. ListPairs:
  133.   PRINT 
  134.   FOR j% = 1 TO NumPairs
  135.     PRINT "Pair";j%,
  136.       Length = LEN(STR$(INT(X(j%))))
  137.       PrintNum$ = RIGHT$(Round$,Length + Places)
  138.       IF X(j%) < 0 THEN PRINT  "-";
  139.       PRINT USING PrintNum$; ABS((X(j%))),
  140.         PRINT  " ",
  141.         Length = LEN(STR$(INT(Y(j%))))
  142.         PrintNum$ = RIGHT$(Round$,Length + Places)
  143.         IF Y(j%) < 0 THEN PRINT  "-";
  144.         PRINT USING PrintNum$; ABS((Y(j%))) 
  145.   NEXT
  146.   PRINT 
  147.   PRINT  "<D>elete a pair  <A>dd a pair  <R>ecalculate  <M>enu"
  148.     
  149.     UpDate:
  150.       UD$ = INKEY$
  151.       IF UD$ = "" THEN UpDate
  152.       UD$ = UCASE$(UD$)
  153.       IF UD$ <> "D" AND UD$ <> "A" AND UD$ <> "R" AND UD$ <> "M" THEN UpDate
  154.       IF UD$ = "M" THEN GOTO PrintMenu
  155.       IF UD$ = "R" THEN GOTO SetDegree
  156.       IF UD$ = "A" THEN AddPair
  157.       
  158.     DeltePair:
  159.       INPUT "Delete which pair"; Which
  160.       IF Which < 1 OR Which > NumPairs GOTO DeletePair
  161.       SWAP X(Which), X(NumPairs)
  162.       SWAP Y(Which), Y(NumPairs)
  163.       NumPairs = NumPairs - 1
  164.       GOTO ListPairs
  165.       
  166.     AddPair:
  167.       NumPairs = NumPairs + 1
  168.       INPUT "X ="; X(NumPairs)
  169.       INPUT "Y ="; Y(NumPairs)
  170.       GOTO ListPairs
  171.        
  172.       
  173.       
  174. PrintPlus: 'print sign before number
  175.   IF POS(1)>60 THEN PRINT 
  176.   IF Flag = 1 AND v(j%)> = 0 THEN PRINT "  +"; 
  177.   IF v(j%) < 0 THEN PRINT  "  -";
  178.   RETURN
  179.  
  180. EndIt:
  181.   WINDOW CLOSE 2
  182.   SCREEN CLOSE 2
  183.   LIST
  184.   END
  185.   
  186. Calculate:
  187.   d = WhatDeg: n = d + 1: d2 = 2 * d
  188.   FOR j% = 1 TO d2
  189.     P(j%) = 0
  190.       FOR k% = 1 TO NumPairs
  191.         P(j%) = P(j%) + X(k%) ^ j%
  192.       NEXT
  193.   NEXT
  194.   P(0) = NumPairs
  195.   r(1) = 0
  196.   FOR j% = 1 TO NumPairs
  197.     r(1) = r(1) + Y(j%)
  198.   NEXT
  199.   IF n = 1 THEN GOTO Jump1
  200.   FOR j% = 2 TO n
  201.     r(j%) = 0
  202.       FOR k% = 1 TO NumPairs
  203.         r(j%) = r(j%) + Y(k%) * X(k%) ^ (j%-1)
  204.       NEXT
  205.   NEXT
  206.   Jump1:
  207.   FOR j% = 1 TO n
  208.  
  209.     FOR k% = 1 TO n
  210.      z(j%,k%) = P(j%+k%-2)
  211.     NEXT
  212.   NEXT
  213.   GOSUB Calculate2
  214.   PRINT : PRINT "degree = "; d
  215.   PRINT "y ="; : Flag = 0
  216.   FOR j% = n TO 1 STEP-1:  IF v(j%)=0 THEN Jump3
  217.     IF j% = 1 THEN 
  218.       GOSUB PrintPlus
  219.       Length = LEN(STR$(INT(v(j%))))
  220.       PrintNum$ = RIGHT$(Round$,Length + Places)
  221.       PRINT USING PrintNum$; ABS((v(j%))) 
  222.       GOTO Jump2
  223.     END IF
  224.     IF j% = 2 THEN
  225.       GOSUB PrintPlus
  226.       Length = LEN(STR$(INT(v(j%))))
  227.       PrintNum$ = RIGHT$(Round$,Length + Places)
  228.       PRINT USING PrintNum$; ABS((v(j%)));
  229.       COLOR 3,0
  230.       PRINT  "x ";
  231.       COLOR 1,0
  232.       GOTO Jump2
  233.     END IF
  234.     GOSUB PrintPlus
  235.     Length = LEN(STR$(INT(v(j%))))
  236.     PrintNum$ = RIGHT$(Round$,Length + Places)
  237.     PRINT USING PrintNum$; ABS((v(j%))); 
  238.     COLOR 3,0
  239.     PRINT  "x^"; RIGHT$(STR$(j%-1),1);
  240.     COLOR 1,0
  241.     Jump2:
  242.     Flag = 1
  243.   Jump3:
  244.   IF POS(1)>60 THEN PRINT 
  245.   NEXT
  246.   q = 0
  247.   FOR j% = 1 TO NumPairs
  248.     q = q+Y(j%)
  249.   NEXT
  250.   m = q/NumPairs: t = 0: g = 0
  251.   FOR j% = 1 TO NumPairs
  252.   q = 0
  253.     FOR k% = 1 TO n
  254.      q = q + v(k%) * X(j%) ^ (k%-1)
  255.     NEXT
  256.   t = t + (Y(j%) - q) ^ 2
  257.   g = g + (Y(j%) - m) ^ 2
  258.   NEXT
  259.   IF g = 0 THEN t = 1: GOTO PrintFit
  260.   t = 1 - t/g
  261.   PrintFit:
  262.     PRINT 
  263.     PRINT INT(t * 100); "% fit"
  264.     IF MenuFlag = 0 THEN HighestDeg = HighestDeg + 1
  265.   RETURN
  266.  
  267. Calculate2:
  268.   IF n = 1 THEN v(1) = r(1) / z(1,1): RETURN
  269.   FOR k% = 1 TO n-1
  270.     a% = k% + 1
  271.     b = k%
  272.     Skip1:
  273.     IF ABS(z(a%,k%)) > ABS(z(b,k%)) THEN b = a%
  274.     IF a% < n THEN a% = a% + 1: GOTO Skip1
  275.     IF b = k% THEN GOTO Skip2
  276.       FOR j% = k% TO n: q = z(k%,j%): z(k%,j%) = z(b,j%)
  277.         z(b,j%) = q
  278.       NEXT
  279.     q = r(k%): r(k%) = r(b): r(b) = q
  280.  
  281.     Skip2:
  282.     a% = k% + 1
  283.  
  284.     Skip3:
  285.     q = z(a%,k%) / z(k%,k%): z(a%,k%) = 0
  286.       FOR j% = k% + 1 TO n: z(a%,j%) = z(a%,j%) - q * z(k%,j%): NEXT
  287.     r(a%) = r(a%) - q * r(k%): IF a% < n THEN a% = a% + 1: GOTO Skip3
  288.   NEXT
  289.   v(n) = r(n) / z(n,n)
  290.     FOR a% = n - 1 TO 1 STEP -1
  291.     q = 0
  292.       FOR j% = a% + 1 TO n: q = q + z(a%,j%) * v(j%)
  293.       v(a%) = (r(a%) - q) / z(a%,a%)
  294.       NEXT
  295.     NEXT
  296.   RETURN
  297.  
  298.  
  299. GraphIt:
  300.   Max = X(NumPairs)
  301.   Min = X(1)
  302.   PRINT 
  303.   PRINT  "Your low X point was "; Min
  304.   PRINT  "Your high X point was "; Max
  305.   PRINT 
  306.   PRINT  "<1> Use these points to graph  <2> Use other points to graph"
  307.   PRINT 
  308.  
  309.   UseWhich:
  310.     UW$ = INKEY$
  311.     IF UW$ = "" THEN UseWhich
  312.     IF VAL(UW$) < 1 OR VAL(UW$) > 2 THEN UseWhich
  313.     IF VAL(UW$) = 1 THEN
  314.       LowX = Min
  315.       HighX = Max
  316.       GOTO FindDiff
  317.     END IF
  318.   
  319.   INPUT "Input the low x value then press RETURN"; LowX
  320.   INPUT "Input the high x value then press RETURN"; HighX
  321.   PRINT 
  322.   
  323.   FindDiff:
  324.   Max = -1E+30
  325.   Min = 1E+30
  326.   Diff1 = ABS(HighX - LowX)
  327.   Diff = Diff1/500
  328.   WhatDeg = 0
  329.   PRINT  "What degree equation do you want to use?";
  330.   FindDegree:
  331.     WhatDeg$ = INKEY$
  332.     IF WhatDeg$ = "" THEN FindDegree
  333.     WhatDeg = VAL(WhatDeg$)
  334.     IF WhatDeg < 1 THEN
  335.       PRINT  "  Invalid entry - try again"
  336.       PRINT  
  337.       PRINT  "What degree equation do you want to use?";
  338.       GOTO FindDegree
  339.     END IF
  340.     PRINT  WhatDeg
  341.     
  342.   GOSUB Calculate
  343.   PRINT 
  344.   PRINT  "Scaling....."
  345.     HighEnd = HighX + (Diff/2)
  346.     FOR Look = LowX TO HighEnd STEP Diff 'find max and min
  347.       YValue = 0   
  348.       FOR k% = 1 TO WhatDeg + 1
  349.         YValue = YValue + v(k%) * Look ^ (k%-1)
  350.       NEXT
  351.       IF YValue < Min THEN Min = YValue
  352.       IF YValue > Max THEN Max = YValue
  353.     NEXT
  354.   CLS
  355.   GOSUB Makegrid
  356.   Scaler2 = 152/(Max - Min)
  357.   Scaler1 = 20 - (Min * Scaler2)
  358.   XCount = 119
  359.     FOR Plot = LowX TO HighX STEP Diff 'plot graph
  360.       YValue = 0   
  361.       FOR k% = 1 TO WhatDeg + 1
  362.         YValue = YValue + v(k%) * Plot ^ (k%-1)
  363.       NEXT
  364.       XCount = XCount + 1
  365.  
  366.       YPlot = (YValue * Scaler2) + Scaler1
  367.       YPlot = 192 - YPlot ' invert
  368.       PSET(XCount,YPlot)
  369.     NEXT
  370.     
  371.     Diff = 500/ABS(LowX - HighX) 'plot the user's points
  372.     FOR PP = 1 TO NumPairs
  373.       NewX = (ABS(X(PP) - LowX) * Diff) + 120
  374.       NewY = 192 - ((Y(PP) * Scaler2) + Scaler1)
  375.       CIRCLE (NewX,NewY),5,5
  376.       CIRCLE (NewX,NewY),6,5
  377.     NEXT
  378.       
  379.    LOCATE 24,23
  380.    COLOR 7,0
  381.    PRINT  "PRESS ANY KEY TO RETURN TO THE MENU";
  382.    COLOR 1,0
  383.    WaitForKey:
  384.      IF INKEY$ = "" THEN WaitForKey
  385.      CLS
  386.      GOTO PrintMenu
  387.  
  388. Makegrid: 'draw garph grid on screen
  389.   FOR Grid = 120 TO 630 STEP 25
  390.     COLOR 2,0
  391.     LINE (Grid,20) - (Grid,172)
  392.     IF (Grid - 120) MOD 100 = 0 THEN
  393.       COLOR 3,0
  394.       LINE (Grid,8) - (Grid,19)
  395.     END IF
  396.   NEXT
  397.   COLOR 2,0
  398.   FOR Grid = 20 TO 172 STEP 8
  399.     LINE (0,Grid) - (620,Grid)
  400.   NEXT  
  401.         
  402.   COLOR 3,0
  403.   FOR Grid = 20 TO 172 STEP 8
  404.     LINE (0,Grid) - (120,Grid)
  405.   NEXT  
  406.  
  407.   COLOR 1,0 'print Y scale
  408.   LOCATE 3,1
  409.   MinMax = (Max - Min)/19
  410.   PrintY = Max
  411.   FOR Grid = 1 TO 20
  412.     PRINT  PrintY
  413.     PrintY = PrintY - MinMax
  414.   NEXT
  415.   
  416.   COLOR 1,0 'print X scale
  417.     Up$ = ""
  418.     WIDTH 80
  419.     AddX = 0
  420.     PadFlag = 0
  421.     UDFlag = 0
  422.     XDiff = ABS(LowX - HighX)/10
  423.     AddX = AddX - XDiff
  424.     LOCATE 1,1
  425.     FOR XX = 1 TO 11
  426.       UDFlag = ABS(UDFlag - 1)
  427.       AddX = AddX + XDiff
  428.       XPrint = AddX + LowX
  429.       XPrint$ = LEFT$(STR$(XPrint),7)
  430.       WHILE LEN(XPrint$) < 12
  431.         XPrint$ = " " + XPrint$
  432.       WEND
  433.       IF UDFlag = 1 THEN
  434.         Up$ = Up$ + XPrint$
  435.         IF PadFlag = 0 THEN
  436.           Up$ = Up$ + "  "
  437.           PadFlag = 1
  438.         END IF
  439.       END IF
  440.     NEXT
  441.     LOCATE 2,5
  442.     PRINT Up$
  443.   RETURN
  444.  
  445. GetIt: 'input routine
  446.   COLOR 7,0
  447.   PRINT  "    --->  ";
  448.   Final$ = ""
  449.   GetItOn:
  450.   User$ = INKEY$
  451.   IF User$ = "" THEN GetItOn
  452.   User$ = UCASE$(User$)
  453.   IF User$ = "C" THEN Final$ = "C": GOTO EndGet
  454.   IF User$ = "D" THEN Final$ = "D": GOTO EndGet
  455.   IF User$ = CHR$(8) OR User$ = CHR$(127) THEN  'backspace
  456.     PRINT User$;
  457.     Final$ = LEFT$(Final$,LEN(Final$)-1)
  458.     GOTO GetItOn
  459.   END IF
  460.   IF User$ = "," THEN GetItOn
  461.   IF User$ = CHR$(13) THEN
  462.     Final = VAL(Final$)
  463.     PRINT 
  464.     IF Final$ = "" THEN
  465.       BEEP
  466.       PRINT  "RETURN with no input - try again"
  467.       PRINT 
  468.       GOTO GetItOn
  469.     END IF
  470.     GOTO EndGet
  471.   END IF
  472.   Final$ = Final$ + User$
  473.   PRINT  User$;
  474.   GOTO GetItOn
  475.     EndGet:
  476.       COLOR 1,0
  477.       RETURN
  478.  
  479. DATA "Least Squares"
  480. DATA " "
  481. DATA "C 1987 to George Trepal   2650 Alturas Rd   Bartow  FL  33830"
  482. DATA "a shareware program - OK to give free but not to sell"
  483. DATA " "
  484. DATA "  This program uses the Least Squares technique to find equations"
  485. DATA "from points.  For example if you can microwave two muffins in"
  486. DATA "25 seconds, three muffins in 50 seconds, and four muffins in 120"
  487. DATA "seconds how long will five muffins take?  To solve the mystery"
  488. DATA "feed the program the point pairs"
  489. DATA " "
  490. DATA "Muffins   Seconds"
  491. DATA "2           25"
  492. DATA "3           50"
  493. DATA "4          120"
  494. DATA " "
  495. DATA "then tell it to calculate.  It'll produce equations of different"
  496. DATA "degrees and judge their quality.  Pick the one of highest quality"
  497. DATA "and feed it 5 to find how long five muffins would take."
  498. DATA "   Of course there are other uses.  You can use it for trend"
  499. DATA "analysis or for graphics programs to find a complex curve that"
  500. DATA "exactly fits your needs.  And it's great to check your homework."
  501. DATA "---> PRESS ANY KEY TO READ MORE <---"
  502. DATA " "
  503. DATA "WARNING!  Equations above the second degree are sometimes tricky."
  504. DATA "Instead of a smooth curve that fits your points they may be snakes"
  505. DATA "that merely intersect them.  Always graph an equation to see what"
  506. DATA "it's realy like."
  507. DATA " "
  508. DATA "Predicted points far away from known points tend to be wrong so"
  509. DATA "use common sense about what you trust.
  510. DATA " "
  511. DATA "If this program can't give you a nice equation from your points don't"
  512. DATA "give up.  There's a lot of math, such as trig functions, this program"
  513. DATA "can't handle.  What It can handle it handles very well but it's not a"
  514. DATA "universal curve finder."
  515. DATA " "
  516. DATA " "
  517. DATA "This is version 1.0"
  518. DATA " "
  519. DATA "Released to public domain September 1987"
  520. DATA " "
  521. DATA "If you like this program please send a contribution."
  522. DATA " "
  523. DATA " "
  524. DATA "PRESS ANY KEY TO CONTINUE"
  525.         
  526.